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Abstract 

In this paper we present a numerical method for hydrodynamic models that arise 
from time dependent density functional theories of freezing. The models take 
the form of compressible Navier-Stokes equations whose pressure is determined 
by the variational derivative of a free energy, which is a functional of the density 
field. We present unconditionally energy stable and mass conserving implicit 
finite difference methods for the models. The methods are based on a convex 
splitting of the free energy and that ensures that a discrete energy is non¬ 
increasing for any choice of time and space step. The methods are applicable to a 
large class of models, including both local and non-local free energy functionals. 
The theoretical basis for the numerical method is presented in a general context. 
The method is applied to problems using two specific free energy functionals: 
one local and one non-local functional. A nonlinear multigrid method is used to 
solve the numerical method, which is nonlinear at the implicit time step. The 
non-local functional, which is a convolution operator, is approximated using the 
Discrete Fourier Transform. Numerical simulations that confirm the stability 
and accuracy of the numerical method are presented. 

Keywords: Classical Density Functional Theory, Phase Field Crystal, 
Compressible Navier-Stokes, Convex Splitting, Finite Difference Methods, 
Energy stability 


1. Introduction 

Solid liquid phase transitions are of great scientific interest. The equilibrium 
properties of this physical process are fairly well understood in the context of 
classical density functional theory (CDFT) of freezing [HE]. This theory char¬ 
acterizes the equilibrium state of a pairwise interacting set of particles in terms 
of the one particle density field. The density function p represents the spatial 
distribution of particles, i.e., the probability of finding a particle at some point 
in space. This function at equilibrium is represented as a minimizer of a free 
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energy, which in turn is a functional of the density. At equilibrium this func¬ 
tion admits two forms of solutions corresponding to a homogeneous distribution 
known as the liquid phase and the inhomogeneous distribution known as the 
solid phase. The inhomogeneous distribution typically consists of peaks located 
on an ordered Bravais lattice representing the probable locations of atoms. This 
approach is attractive for the reason that the equilibrium solid phase carries the 
information about the lattice symmetries. A non-equilibrium time dependent 
theory that characterizes the non-equilibrium distribution has the potential to 
capture lattice-dependent anisotropic effects such as defect formation and mi¬ 
crostructure evolution. Such models have many potential applications in the 
area of materials modeling. Consequently, the development of time dependent 
models has been the focus of recent research 00000001m. 

Given a free energy functional ^[p] the simplest dynamical equations for the 
evolution of the density field are to use gradient dynamics on the free energy 
surface mm- While this approach is appealing it has several shortcomings 
including the inability to capture the inertial effects and the effect of flow on 
the phase transition. Hydrodynamic coupling has been incorporated in models 
involving colloidal suspensions mm and in models that describe the freezing 
of a dense hard sphere gas m • In all these approaches the hydrodynamic 
model takes the form of a compressible Navier-Stokes-like equations where the 
pressure (p) is defined by the equation of state Vp := pV^. This can be easily 
seen from a thermodynamic point of view through the Gibbs Duhem relation 
dp = pdp where p is the chemical potential taken to be the variational derivative 
of the free energy [ 8 ]. Models for hydrodynamic coupling with solid-liquid phase 
transitions thus take the form 


dtp + V • (pu) = 0 , 

d t (pu) + V • (pu 0 u) = — pV + 7 V 2 u, 


(1) 


where p(r, t) : ft x [0, oo) —>> M d , u(r, t) : Vt x [0, oo) M d , Q C M d , ^[p] is a 
functional of p which may be local or non-local in nature, 7 > 0 is a viscosity 
coefficient and d is the dimension of the system. It is easy to see that the model 
in Eq. 0 admits an energy 

£\pM = \ J q P I u I 2 dr + F[p], (2) 


such that 


f <0. 

dt 


(3) 


The proof of Eq. 0 is presented in |Appendix A| The long time equilibrium 


for the above system corresponds to 5 p J-[p] = 0, u = 0. Thus the equilibrium 
properties of the system are determined by the functional J~[p\ which represents 
the free energy functional. The specific functional form of the free energy J-[p\ 
determines the theory and the phase diagram of the model. One must appeal 
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to CDFT or a similar approximation for the specific form of the free energy. 

It is worth noting that the time dynamics of the system however is much 
richer than a standard gradient descent model m • The first set of simulations 
using this type of model was presented in m- It was observed that while the 
total energy £[p,u] and the free energy J-[p\ were non-increasing functions of 
time, the kinetic energy and the individual components of the free energy were 
not. However no details about the numerical method were given in pm- Here, 
we present the first energy-stable method to solve such coupled systems. 

The physical properties and relevance of the model lies in the free energy 
functional J-[p\ used. A large number of models for the free energy exist in 
literature. These models include CDFT BE] derived in the framework of equi¬ 
librium statistical mechanics, as well as phase field crystal (PFC) models pT] , 
which are derived using phenomenological theories motivated by CDFT. The 
energy in CDFT is typically a non-local functional of the density while the en¬ 
ergy in the PFC is a weakly non-local functional that depends only on gradients 
of the density. In addition, in CDFT the ideal gas part of the free energy con¬ 
tains a logarithm. In CDFT, the combination of the logarithmic part of the 
free energy with the nonlocal interaction energy induces much sharper peaks 
in the density field than the peaks obtained using the PFC model where the 
logarithmic term is replaced by a polynomial and a gradient approximation is 
used for the interaction energy (as described in Appendix B). Thus, the CDFT 
system requires a finer grid to resolve the dynamics. 

These free energy models have been very popular and capture phase tran¬ 
sitions at the atomic length scales, but on diffusion time scales, thus enabling 
the simulation of microstructure evolution at much longer time scales than can 
be captured using molecular dynamics or other stochastic approaches m- This 
makes the hydrodynamic theory considered here an extremely valuable tool in 
investigating non-equilibrium effects such as the effect of flow on the phase tran¬ 
sition. Thus developing an efficient and stable numerical method for this class 
of models is of critical importance. 

In this work we use the convex splitting framework introduced by Eyre [12] 
to develop stable numerical schemes. The convex splitting framework has been 
successfully applied to the PFC-based models pH13 ED3 DU HE and to non¬ 
local CDFT-like models [18j [19] . However these approaches consider conserved 
gradient descent models without hydrodynamics. The current literature on 
hydrodynamic models driven by free energy gradients mainly considers incom¬ 
pressible multiphase fluids using the Cahn-Hilliard free energy (see [20] for a 
comprehensive review). We note that in GB, the incompressible Navier-Stokes 
system was coupled with the PFC model. A similar approach using a phe¬ 
nomenological PFC-like free energy driven hydrodynamic model was presented 
in J22j however no numerical methods were presented. While there has been 
much work in the literature on the development of numerical methods for com¬ 
pressible fluids driven by van der Waals and Cahn-Hilliard-like free energies, to 
the best of our knowledge the development of accurate and energy-stable nu¬ 
merical methods for the case of compressible fluids driven by gradients of PFC 
and CDFT free energies, which are higher order and nonlocal (CDFT), has not 
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been addressed. 

The paper is structured as follows. In Section [2] we present the description 
of the model and the associated free energies. In Section [3] we present the basic 
ideas of convex splitting and outline the need for a comprehensive treatment for 
the various free energy models discussed in this work. This section also presents 
the basic energy estimates required to develop an unconditionally energy stable 
method for the hydrodynamic theory of freezing. An unconditionally energy 
stable conservative discretization of time is presented in Section [4| where the 
properties of scheme are proven. The corresponding fully discrete numerical 
method is presented in Section [5j Finally the numerical results are presented in 
Section [6j where the properties of the numerical method are explored numeri¬ 
cally. In Section [7| conclusions are drawn and future work is discussed. Details 
of the model derivation, definitions of the discrete operators, and the nonlinear 
multigrid method used in the simulations are presented in the Appendices. 


2. Hydrodynamic Model and the Associated Free Energies 


In the reminder of the paper we will restrict ourselves to 2 dimensions. How¬ 
ever we note that the ideas, numerical methods, theorems and proofs extend to 
3 dimensions in a straightforward manner. We will consider a spatial domain 
D := [0, L x ) x [0, Ly) with periodic boundary conditions. The results however ex¬ 
tend to other boundary conditions such as homogenous Dirichlet (no-slip) solid 


wall boundary conditions and Neumann boundary conditions (see Appendix 
E.2). The hydrodynamic equations Eq. 0 are first rewritten in the primitive 


variable form for convenience as follows: 


dtp + V • (pu) = 0, (4) 

p(d t u + u- Vu) = -pV + 7 V-P. (5) 

The dissipation tensor V is defined as 

D ;= (Vu + Vu T ) - (V • u)I, (6) 

where I is the identity tensor of rank 2. In this work we will consider the 
following two specific models for the free energy. 


2.1. Classical Density Functional Theory 

Classical density functional theory typically relies on a free energy that is a 
non-local functional of the density p. A typical form for the free energy is given 


by 

Fcdft[p\ = Fid[p\ + Fex [p \, 

(7) 

where 

Fid[p\ := [ p(ln(p) - l)dr, 

Jn 

(8) 
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is the ideal gas part of the free energy (e.g., describes mixing) and 

Fx[p] ■= J^P( J * P ) dr > (9) 

which is the excess free energy functional (e.g., describes interactions), with 

( J*p)~ [ J(x-y)p(y)dy. (10) 

Jn 

The convolution kernel J represents different quantities in different approxi¬ 
mations of the density functional theory. For example, in the context of the 
local density approximation J(r — r') = J(\ r — r' |), represents the pairwise 
inter-atomic potential energy associated with two particles that located at r and 
r' respectively (see HD- In the context of the Ramakrishnan-Yousseff density 
functional theory the quantity J represents the two particle direct correlation 
function [23] . While the convolution kernel has different physical meanings for 
different physical approximations the form of the excess free energy is consis¬ 
tent across a large class of density functional theories and is worth addressing 
in general. 

For the reminder of this paper we will assume that the convolution kernel J 
is a sufficiently regular and fl-periodic function such that 

• J = J c — J e , where J c , J e are sufficiently regular, ^-periodic and pointwise 
non-negative. 

• J c and J e are even, i.e., J a (—r) = J a ( r) where r £ M 2 , a = c, e. 

• In J ( r )dr < 0. 

It is worth noting that the above assumptions are in general true. For example, 
the nontrivial assumption that the kernel is even follows from radial symmetry 
of the interaction potentials for both the local density approximation and the 
Ramakrishnan-Yousseff formalism [TJ. 

This free energy satisfies a convex splitting of the form Tcdft •= Fcdft, c — 
FcDFT,e where Tcdft,c and TcDFT,e are convex functionals defined as: 

FcdftApI '■= [ {yo(ln/O — 1) + (Je * l)p 2 } dr; 

Jn (n) 

FcdftAp\ ■= J R/°(J * P) + (Je * 1 )p 2 \ dr. 

Here we say a functional .F[p] is convex in a Hilbert space H (consisting of 
sufficiently regular functions) if for any p G H that lim e ^o £2 T\p + ev] > 0 for 
all v £ H . This convex splitting is not unique and other forms of convex splitting 
can be obtained (see HU US]). This specific convex splitting is similar to the 
one proposed in m as it has certain advantages, which will become apparent 
in the later sections. It is worth noting that in m the authors considered the 
case where (J* 1) > 0. However in CDFT (J* 1) is usually negative as assumed 
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here. For example in the Ramakrishnan-Youssef [23] CDFT model (J * 1) is 
related to the isothermal compressibility of the liquid and is negative pQ. 


2.2. Phase Field Crystal Model 

The phase field crystal (PFC) approximation was first proposed by Elder et 
al [24]. While the PFC was first derived using a phenomenological approach, 
the PFC has since been re-interpreted as an approximation to dynamic density 
functional theory (DDFT), which is a time-dependent counterpart to CDFT 
[Him m nu . Typically, PFC formulations use a scaled version of the deviation 
of the density from a reference density rather than the density field itself (see 
Appendix B). Asa result, the density is not always non-negative in PFC models. 


Here, we use a formulation of the PFC that uses the particle density field directly 
and maintains its non-negativity. 

The PFC energy functional we consider is given by 


7ffc(i,) =L {n (" “ I) + ?(' , -0 


TT ( P-o ) -|Vp|' + j(A„)' *, (12) 


with 


SFpfc 
Sp 


= o\P 


p- 2 ) + 2A P- 


(13) 


where a is a free parameter that determines the strength of interactions between 
the particles and hence the elastic constants of the solid phase. We refer the 
reader to | Appendix B for details regarding the relation of the above free energy 
and more standard formulations of the PFC energy. We note that this free 
energy admits a convex splitting of the form Tpfc ~ ?pfc,c ~ ^PFC,ei where 
Fpfc,c and PpFC,e are convex functional defined as: 


and 



(14) 


(15) 


3. General Framework of the Convex Splitting Schemes 

In this section we discuss in considerable generality the ideas involved in 
convex splitting schemes for the class of models described so far. Previously, 
Wise et al in m considered a class of local free energies like the PFC model that 
satisfies the form T := T c —T e ^ where F e [p\ and T c [p\ are convex functionals of p 
with a specific functional form. It was assumed that T e [p\ and T c [p] are integrals 
of energy densities / c (p, d x p, d y p, Ap) and / e (p, d x p, d y p, Ap) respectively. It 
was further assumed in m that the functions f c and f e are convex in their 
arguments respectively. The authors derived an energy estimate for this special 
case as summarized in the following theorem. 
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Theorem 3.1. Given a free energy F[p] that admits a convex splitting of the 
form F := F c — F e where F e [p\ and F c [p\ are convex functionals of p of the 


form 

r c [p\ = 

/ fc(p,d x p,d y p,Ap) 

Jn 

(16) 

and 

TM = 

/ fe{p,d x p,d y p,Ap) 

Jn 

(17) 


such that f c and f e are convex in their arguments. Further if p is sufficiently 
regular and periodic in Cl and d x p , d y p are also periodic in Cl, then 


F[(p] - J U [lp] < (Safety] - Safety), <j> - Ip ) 2 , (18) 


where (•, -)2 is the usual L 2 inner product and 5 P represents the variational 
derivative with respect to p. 


Remark 3.2. Note that the estimate Eq. (18) can be immediately translated 
into an energy estimate. Consider a discretization of time t k G [0, T] such that 
t k := 0 + ks where s > 0, s G M and k G Z + . Then we have the energy decay 
estimate 


F[p k+1 ] - T[p k ] < (6 P F c [p k+1 ] - 5 p F e [p k ],p k+l - p k ) 2 , (19) 

where p k := p(x, t k ). This energy decay estimate can be used to show energy 
stability for convex splitting schemes. 


The energy estimate however is limited to the class of free energy functionals 
considered in Theorem 3.1 The CDFT free energy presented in Eq. (|7|), which is 
a nonlocal functional of the density field, is an example of a free energy functional 
out of the scope of Theorem |3.1| Nonlocal free energies were considered in 
mm, where a convex splitting scheme was presented taking advantage of the 
specific form of the CDFT free energy functional. 

In what follows we consider a general approach to convex splitting schemes 
and the associated energy decay estimates that seamlessly include a much larger 
class of free energies including the PFC and CDFT. Before we proceed further 
we define the notion of a proper convex splitting: 


Definition 3.3. Consider a functional F[p\ : H C L^ er {Cl) M, where H C 
Lp er (Cl ) is a Hilbert space of sufficiently regular periodic functions p : Cl — >• M 
with Cl = [0, Lx) x [0 ,L y ). F[p] is said to admit a proper convex splitting if 
T := T c — T e such that 


de 2 


F a [p + ev\ > 0 


Vp, v G H,e Gl, a G {c, ej. 


( 20 ) 
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Note that the definition above is stronger than a simple convex splitting, where 
T c and T e are convex functionals, i.e., 

d 2 

lim ~T 7 ;J = 'a[p + erl > 0 Vp, r G H, e G M, a G ic, e}. (21) 

e-^o de^ 

If T := T c — T e is a proper convex splitting of T then the T c and T e are 
indeed convex functionals. However the converse is not necessarily true. Taking 
advantage of the proper convex splitting of T we can obtain the energy decay 
estimate presented in Theorem |3.1| 

Theorem 3.4. Given a free energy J-[p\ that admits a proper convex splitting 
of the form T := T c — T e , where J~ e [p\ and J~ c [p\ are functionals of p such that 

d 2 

[p + ev] > 0 Vp, v G H, e G M, a G {c, e}, (22) 

de z 

then 

F[4>\ - F[$\ < (b<t)Fc[4>\ - 0 - '0)2, (23) 

where (•, 02 is usual L 2 inner product and S p represents the variational 
derivative with respect to p. 

Proof. The proof is presented in |Appendix C| 

In order to demonstrate the power of this decay estimate let us first consider 
the simple gradient model : 

dtp = V 2 S p T. (24) 

Using the energy estimate in Remark. |3.2| we can construct an unconditionally 
energy stable conservative scheme for gradient dynamics with free energies that 
satisfy a proper convex splitting. This is summarized in the following theorem. 

Theorem 3.5. Suppose that Q = [0 ,L X ) x [0 ,L y ) and p k : U —>> M 2 and u k : 
U —>> M 2 for all k > 0, A; G Z + are periodic and sufficiently regular. Then the 
solution to the scheme 

P k+1 -P k = sV 2 (5 p F c [p k+1 ] - 5 p Fe[p k }), (25) 

conserves mass, i.e., 

(y+\i ) 2 = (y,i) 2 , vfc > 0. (26) 

If in addition T admits a proper convex splitting such that 

F [ 4 >] - - VeM, 0 - 0)2, (27) 

then the scheme is unconditionally energy stable and 

F[p k+1 ] < Hp k ] Vfc > 0, (28) 
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for any s > 0 . 


Proof. Conservation follows by integrating Eq. (25) over Q to obtain 

(p k+i ,i) 2 = ( P k ,i) 2 . 


Choosing <f = p kJrl and if = p k in Eq. (27) we have the estimate : 

F[p k+1 ] - F[p k ] < (<5 P .F C [/ +1 ] - S p r e \p k ],p k+1 - p k h- 


Now using Eq. (25) to replace p k+1 — p k we have 


7[p k+1 ] - T[p k ] < -s || V {5 p T c [p k+1 \ - S p T e [ P k ]) |||< 0 . 
In the above expression || • |||:= (•, -) 2 - 


( 29 ) 


(30) 


(31) 


Remark 3.6. The above theorem seamlessly includes both the local and non¬ 
local free energy functionals discussed in Section [1| i.e., the cases of PFC and 
CDFT, thus serving as a stronger result than the one presented in Wise et al 
m- Further Theorem \S. 1\ follows immediately from Theorem \S.J\ To see this 

and T e satisfy 
The proof of 


note that convexity of f c and f e in their arguments implies T c 
Eq. [20 1). Hence T = T c — F e forms a proper convex splitting. 


Theorem 3A now follows from Theorem\S.J\ and can be treated as a corollary. 


Remark 3.7. D. Eyre Tfflj is credited widely for the convex splitting approach 
which was originally applied to the Cahn-Hilliard equations. This is the first 
instance to the best of our knowledge where the energy decay was characterized 
solely in terms of a convexity property without specific functional forms, which 
enables convex splitting arguments to apply to a wider class of free energies, i.e., 
the ones that satisfy what we call a proper convex splitting. In a similar manner 
one can construct convex splitting schemes for the hydrodynamic models that 
form the subject of this paper. This is presented in the following section. 


4. Semi-discrete Energy Stable discretization of the Hydrodynamic 
Model 

An unconditionally energy stable time discretization of the model in Eq. Q 
takes the form (space remains continuous): 


p-~ ~ P k =-sv-(y u fe+ 5), 

k 

-p k UJ k X U fc+ 5 - ^-V I U fe+1 |2 -" k 't7" k+1 


p k (u k+1 - u fc ) = 


-p k Vp, K 


i rT 1 

,fc+ b p 

- y v ' “ 

+ 7 V • V k+ 5 

M fe+1 = 5 p T c [p k+1 ] - 6 p T e [p% 


(32) 

(33) 

(34) 


where we have used the relation u • Vu = wxu + \ V | u | 2 and w = V x u. In 
the above equations the superscripts ( k and k + 1) stand for the discretization 
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of time as in the case of Remark 13.21 and 

u fc+ 5 := i(u fc+1 +u fc ), (35) 

and 

pfc+l , = (vu fc+ 5 + (Vu* + 5) r ) - V • U fc+ 5I. (36) 

It is easy to see that V • V k+ ^ = V 2 u. Note that this discretization respects 
the proper convex splitting of the free energies, i.e., the convex part of the free 
energy is treated as an implicit term and the concave parts are treated as explicit 
terms. Now we have the following theorem. 

Theorem 4.1. Suppose that £2 = [0 ,L X ) x [0 ,L y ) and p k : £2 M 2 and u k : 
£2 —>> M 2 for all k > 0, fc G Z + are periodic and sufficiently regular. Then the 
solution to the scheme Eq. conserves mass, i.e., 

(/ +1 ,1) 2 = (y,l) 2 , Vfc > 0. (37) 

If in addition T admits a proper convex splitting such that 

F[<i>] - - VO 2 , (38) 

then the scheme is unconditionally energy stable with respect to the total energy 
£[p, u] (defined in Eq. [1|) and 

£[p k+1 ,n k+1 ] < £[p k , u fe ] Vfc > 0, (39) 


for any s > 0. 


Proof. Integrating Eq. (32) over £2 we immediately have mass conservation 
(p fe+1 , l) 2 = (p fe , l) 2 . Choosing <f = p k+1 and if = p k in Eq. (38) we have the 
following estimate on the free energy : 

T[p k+1 \ - T[p k } < (S P T, ' c [p k+1 \ ~ S p T e \p%p k+1 - p k ) 2 . (40) 

This estimate immediately translates to an estimate for the total energy (see 
Eq. ©) as follows : 

£[p k+1 ,u k+1 ]-£[p k ,u k ] = ±(p k+1 u k +\u k+1 ) - i (/u fe ,u fc ) 

+T[p k+1 ] - T[p k ] ( 

< t (p k+1 u k+1 , u k+1 ) - I (p k u k , u k ) 

+(p k+1 - p k ,p k+1 ), 


where we have used Eq. (40) in in step 2 along with the definition of p k+1 (Eq. 
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(34)). Multiplying Eq. (33) by u fc+ 2 and integrating over ft we have 

1 f p k+1 I u fc+1 I 2 dr - l f p k I u fc I 2 dr - J [ (p k+1 - p k ) | u fc+1 | 2 dr 

2 Jn 2 J n 2 J n 


/c+l |2 


= -si [ p k u fe+ 5 • V I u 
2 Jn 

-s / p k u fe+ 5 • vy +1 d: 

Jn 

+s 7 J u fe+ 5 • (v-D fe+ 5)dr. 


Now we note that 


-1 [ (p k+1 - p k ) I u fe+1 I 2 dr = l [ sV • (/u fe+ 5) I u fe+1 ' 2 

2 Jn % Jn 

= -si I p k u fe+ 5 • V I u 
2 Jo 


/c+l 12 


dr, 


dr, 


(42) 


(43) 


where we have used the continuity equation Eq. (32) in the first step and 


integration by parts in the second step. Further we have 

—5 J /u fc +5 • Wp k+1 dr =s j V ■ (p k u k+ i^J p k+1 dr, 

= - [ (p k+1 - p k ) p k+1 dr, 

Jn 

where we have used integration by parts in the first step and the continuity 


(44) 


equation Eq. (32) in the second step. Now combining the relations Eq. (43) 


and Eq. (44) with Eq. (42) and simplifying the expression we get 


i(p fc + 1 u fc+1 ,u fc+1 ) - l (p k u k ,u k ) 


+(p k+1 - p k , p, k+1 ) = SJ 


' [ U fe+ 2 

Jn 


. V-P fc +2dr, 


(45) 


= : V k+ ?dr, 

2 Jn 


where we have used the identity 

* k +h. . V7 . 'Tik+h.rlr — 


[ u fe+ 3 • V-V k+ ^dr = [ V k+ i :V k+ idr, 

Jn 2 J Q 




(46) 


in the second step. This expression is proven in | Appendix D| Finally using Eq. 
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(45) to simplify the right hand side of Eq. (41) we have 


£[p k+1 , u fc+1 ] - £[p\ u fe ] < \{p k+l u fc+1 , u fe+1 ) - t (p k u k , u k ) 

+ (p k+1 - p k ,p k+1 ), 

= [ £> fc+ 3 : P fc+ 5rfr, 

^ J Q 

< 0. 


( 47 ) 


Hence the total energy of the system is non-increasing, regardless of the time 
step s > 0 i.e., £[p k+1 , u fc+1 ] < £[p k ,u k ]. Such schemes are called uncondition¬ 
ally energy stable. 

The above theorem is one of the main results in this paper. The discrete 
time continuous space method is energy stable for the CDFT and the PFC free 
energy models presented in the previous section. It is worth noting that the 
proof of the above theorem relies entirely on the free energy satisfying a proper 
convex splitting. Thus the theorem holds for a large class of continuum free 
energy models where the free energy admits a proper convex splitting. In what 
follows we will present the fully discrete version of the method. 


5. Fully Discrete Numerical method for Hydrodynamic Models 


A fully discrete energy stable discretization of Eq. 0 is written as : 


i +1 - Pij = -sv h ■ (p k j u yh 


pijfaij 


u t)= s 


44 x 4A - ttW i u y 1 1 2 -4w44 


+7A.UJA 


= W4 +1 ]-W4h 

4 = (v„ x 4 .), 


(48) 

(49) 

(50) 

(51) 


where p G C mxn periodic, J 7 , T c , T e : C mX n — > M are fully discrete free en¬ 
ergy functionals such that .F[p] = Tdp] ~ Fe[p\ is a proper convex splitting of 
J-'lp]. The discrete function spaces and the corresponding discrete differential 
operators and norms are defined in | Appendix E| 


5.1. Unconditional Energy Stability and Conservation 

In order to show that the fully discrete scheme presented in Eqs. (48) - (51) 
is unconditionally energy stable and mass conserving we present the following 
general theorem which is independent of the form of the free energy. The result 
however takes advantage of the energy estimate from the proper convex splitting 
to show unconditional stability. 
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Theorem 5.1. 

0 are periodic , 
conserves mass , 


Suppose that p k G Cmxn and u k 
then the solution to the scheme 


^ Cmxn x Cmxn far all k ^ 

Eq. ~ & p k+ \ 


u 


k+1 


i.e.. 


{p k+1 || 1) = {p k || 1) , Vfc > 0. (52) 

If in addition T : = T c ■— J~ e satisfies the inequality 


F[4>\ - - 6i,?e [i/>] H-tP), (53) 


then the scheme is also unconditionally energy stable, with discrete energy 

S[p k ,n k ]= l -{p k n k \\n k ) + T[p k ] . (54) 


Thus we have 

£[p k+ \u k+1 ] < £[p k , u fe ] Vfc > 0, (55) 

for any s > 0. 


Proof. The proof of the theorem is similar to Theorem 4.1 and is presented in 
|Appendix F Note that (• || •) is an appropriately defined discrete inner product 
(see Appendix E). 


Now in order to construct a fully discrete convex splitting scheme that is 
unconditionally energy stable and mass conserving all we need is a fully discrete 
proper convex splitting of the free energy. 


5.2. Proper Convex-Splitting for Fully Discrete CD FT Free Energy 

The fully discrete free energy functional the CDFT model Tcdft • C mxn 
M is given by : 


Tcdft\p\ = (p || ln(p) - 1) - * P)- ( 56 ) 

where J, J c , J e G C mxn are the restrictions of the convolution kernel given by 
Ji,j •= J(xi,Xj). In the above equation the periodic discrete convolution is 
defined as: 

m n 

{J * p)i,j •= h 2 ^ ^ Jk,lPi-k.j-h (57) 

k l 

where h is the uniform grid spacing for Cmxn (see 
We now have the following energy estimate. 

Lemma 5.2. Suppose p G Cmxn is periodic then the energies 

FcdftM : = (p II l n (p) - 1) + {Je * 1) I P ||i 5 (58) 

and 

J~CDFT,e [p\ := {Je * 1) || P ||| {p || J * p), (59) 


Appendix E for definitions). 
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are convex in p and Tcdft •= Tcdft,c ~ J~CDFT,e is a proper convex splitting 
of the Tcdft • The gradients of the energies are 


dpT cdft,c[p\ = ln(p) + 2(J e * 1 )p, (60) 

dpTcDFT,e[p\ = 2(J e * l)p + J * p, (61) 

and the energy satisfies 

m - Fbl>] < {SpFc[4>\ - SpTebl>] II 0 - # (62) 


Proof. It is easy to verify that 


-^FcdftAp F q4>\> 0 


and 


-j^F C DFT,e[p F q4>}> 0 


for any periodic </> G C mX n 5 ^ £= M. The inequality in Eq. (62) follows from the 
proper convex splitting of Tcdft • 

Corollary 5.3. Suppose that p k G Cfhxn and u k G Cfhxn x Cfhxn for all k > 0 
are periodic and the discrete free energy is defined by Eqs. 
scheme in Eqs. 
for s > 0. 


51|) zs mass conserving and unconditionally energy stable 


Proof. The proof follows from Theorem 5.1 and Lemma 5.2 


5.3. Proper Convex-Splitting Scheme for Fully Discrete PFC Free Energy 

The fully discrete free energy functional the PFC model Tpfc : C mxn —> M 
is given by : 


r P Fc[p\ L || p _ 3 || 4 + :| || p „ 3 ||a + 1 || || 2 + (p || A ,p), (63) 


where the last term is introduced to approximate — || Vp || 2 in Eq. (12) and 
the discrete norm || • ||4 is defined in Appendix E The corresponding energy 


estimate is now given by the following lemma. 


Lemma 5.4. Suppose p G C mxn is periodic and A^p is also periodic then the 
energies 

TpfcAp] ■= ^ II P- § Ik +f || p- \ || 2 II A h p || 2 , (64) 

and 

FpfcApI '■= ~ (P II A hp ), (65) 
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are convex in p and Tppc J~pfc,c ~ J~PFC,e is a convex splitting of the 
Tpfc • The gradients of the energies are 

5pTpFC,c[P\ = +a ( P ~'^) + A ^’ ( 66 ) 

bpFpFCT,e[p\ = —2A hp, (67) 

and the energy satisfies 

m - Tty] < (SpTM - 5 p E e ty] || 0 - VO- (68) 


Proof. The proof is similar to that of Lemma [5T2] and we omit the details for 
brevity. 


Corollary 5.5. Suppose that p k G Cfhxn an d u k E Cfhxn x Cfnxn f or all k > 0 
are periodic and the discrete free energy is defined by Eqs. (63)- (65), then the 
scheme in Eqs. f7#| ) - [5l\ ) is mass conserving and unconditionally energy stable 
for s > 0. 


Proof. The proof follows from Theorem 5.1 and Lemma 5.4 


6. Numerical Results 


To solve the nonlinear system at the implicit time level, we use a standard 
non-linear Full Approximation Scheme (FAS) (see [27]). The implementation of 
the prolongation and restriction operations for the multigrid solver are presented 
in detail in |28] in the context of a linear problem. The Gauss-Seidel smoothing 
scheme used in the FAS method is presented in | Appendix G.l and |Appendix| 
|G.2| respectively for the DFT and PFC models. Further we note that we use 
one smoothing iteration per level in all the simulations presented in this paper. 
Because the density field takes the form of sharp Gaussian like peaks on a lattice 
with near vacuum in the interstitial spaces as observed in m in the case of the 
CDFT models, two modifications to the method are introduced: 


1. The first is to introduce a projection step for the density field. We use 
a projection step where we replace p by its positive projection yj p 2 + S 2 . 
Here S is a regularization parameter chosen to be 10“ 10 in the case of the 
CDFT model. This allows us to handle the vacuum seamlessly and ensure 
that the log(p) term does not receive a negative input in the vacuum 
region. The parameter S is set to 0 in the PFC case. 


2. The second modification is to damp the Gauss-Seidel iterations and slow¬ 
ing down the approach to the solution (see Appendix G.3). This pre¬ 
vents the method from overshooting the vacuum and producing negative 
densities. Although preliminary simulations indicate that damping alone 
might be sufficient to maintain positivity of the density the combination 
of damping and projection allows us to use minimal damping and faster 
convergence. 


15 













In the reminder of this section we present numerical evidence for convergence 
and energy-stability of the schemes. We also present some simulations illustrat¬ 
ing the mechanisms of the solid-liquid phase transition in the model and the 
ability of driven flows to steer the system away from equilibrium. We use the 
following model convolution kernel: 

J(x) = a /2 exp - ex P (-^ 2 x2 ) - ( 69 ) 

where v is a constant spatial scaling. The specific form is chosen to model a 
two particle direct correlation function [lj with a single peak in Fourier space. 
The parameter v is chosen to recover a phase transition qualitatively similar to 
the one presented in m- The function J e is chosen to be the function : 


f — J(x) if J(x) < 0 
| 0 otherwise 


(70) 


J c is now defined by the relation J c = J + J e . All simulations presented in this 
work use a residual tolerance of 10“ 12 for the PFC model and 10 -14 for the 
CDFT model. The viscosity coefficient 7 is set to 7 = 2 in all simulations. 


6.1. Convergence in time 

In order to estimate the convergence rate we evolve the system using same 
initial data with increasingly finer grid resolutions. The time step refinement 
path is chosen to be s = 0.025ft 2 . This choice is not motivated by a time step 
restriction as the scheme is unconditionally energy stable. Since the scheme is 
first order accurate in time and second order accurate in space the refinement 
path of s = 0.025ft 2 should predict a global error of order 0(h 2 ). This is verified 
numerically below. 


6.1.1. Hydrodynamic CDFT Model Convergence Test 

For the case of the hydrodynamic CDFT model we use smooth initial con¬ 
ditions analogous to those previously used for the PFC models in mm - The 
density field used is 


p(x,y,t = 0 ) = 


_ o/7r(x + 10)\ o{7r(y + 3) 
+ 0.02 cos 2 ( —- ] cos 2 1 J 


32 


32 


nm • 2 (^x\ . 2 ( 4tt(j/-6)^ 

■°° ism ("32" ) s,n ( 32 ) ’ 


(H) 


and p 0 = | x 0.6 which is the density corresponding to packing fraction 0.6 
(see uni). The initial velocity field is set to u = 0, i.e., a stationary field. The 
parameter v = 1 in this test case. 

To estimate the convergence rate we start with a grid spacing h and time 
evolve the system to tf = 10 using three different grid spacings ft/ 2 , ft and 2 ft. 
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Given solutions p h on a grid of spacing /i, we define the Cauchy error between 
two solutions p h and p h / 2 as 


h]h /2 
e ij 



hi 2 

P2i,2j- 


The convergence rate is then defined as 


log 2 


VII e ^ /2 IbV 


(72) 


(73) 


The results for grid sizes 16 2 , 32 2 , 64 2 ,128 2 , 256 2 and 512 2 are characterized in 
Tables [l][3] for the density field and the two components of the velocity field. The 
results suggest that the method is in fact first order accurate in time and second 
order accurate in space. The energy evolution shown in Figure [l] demonstrates 
that the total energy of the system is in fact non-increasing. The figure also 
shows a striking feature of the model observed in uni that the model does not 
preserve monotonicity of the kinetic energy or the free energy of the system. 


6.1.2. Hydrodynamic PFC Model Convergence Test 

We use a smooth initial condition analogous to the one previously used for 
the CD FT convergence test given by 


p(x, y,t = 0) = 


V 3 


—0.07 — 0.02 cos 


'2n{x — 12 ) 


+0.02 cos 2 


' 7 t(x + 10 ) N 
32 


32 


sm 


r(2/ + 3) 


32 


-0.01 sin 


2 [ s i„ 2 ( 4lrfa ~ 61 


32 J \ 


32 


"2 ?r (y - 1) \ 

V 32 J 

) 

1 

2 5 


(74) 


where ( x,y) T [0,32) x [0,32). This corresponds to the initial conditions 

used in the standard PFC model convergence test case presented in ii. 
With this as the density field and stationary initial velocity field we choose 
parameter e = 0.025 and time evolve the system to tf — 10. These parameters 
predict a non-uniform steady state corresponding to a solid phase (demonstrated 
in Section 6.2). The other parameters are chosen to be identical to the CDFT 


model presented in the previous section. The results summarized in Tables [4][6] 
suggest that the method is first order accurate in time and second order accurate 
in space for the PFC free energy as well. The total energy is also observed to be 
non-increasing (see Figure [2|. It can be seen in Figure [ 2 ] that the hydrodynamic 
PFC model also does not preserve the monotonicity of the free energy or the 
kinetic energy. 


6.2. Simulation of Solid Liquid Phase Transition 

Next we simulate freezing of a super-cooled liquid using the two models. In 
both cases we start the simulation using a random perturbation of a constant 
density field p = p{ 1 + O.I 77 ), where p is a random number in [0,1]. The average 
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density p is chosen to lie in the solid phase. The simulation of the CDFT model 
is performed with p = | x 0.6 corresponding to packing fraction 0.6. The 
scaling parameter is chosen to be v = 2.362. The simulation is performed on 
a domain Cl = [0,17) x [0,17) with time step s = 0.02 with 128 x 128 uniform 
discretization. The discrete Fourier image J(k x ,k y ) is shown in Figure [ 3 J where 

)e x p (75) 

£=1 j=1 ' ' 

and (, X£,yj ) = (£h,jh) and (. k x ,k y ) belong in a 128 x 128 discretization of 
[— 7 r, 7 r) x [— 7 r, 7 r). Observe that there exists wave numbers k = ( k x ,k y ) T such 
that 1 — pJ( k) < 0 so that the homogeneous density p is unstable, as described 
in m- The time evolution of the density and velocity fields are shown in Figure 
[4] The figure shows how the flow field drives the density field to form the peaks 
located on a lattice that represent the solid phase. The corresponding energy 
evolution is shown in Figure |5j It is seen that the total energy is a non-increasing 
function. Figure [5] shows that during the liquid to solid phase transition the 
ideal gas part of the free energy increases because the system is acquiring order. 
Concomitantly, the excess part of the free energy decreases. The total Helmholtz 
free energy and the total energy of the system are observed to be decreasing 
functions of time. Once again the kinetic energy is a non-monotone function of 
time showing a sharp increase in kinetic energy at approximately time t ~ 475, 
which corresponds to the sharp decrease in the free energy and total energy 
of the system. This corresponds to the fast relaxation and pattern formation 
during the freezing process (see Figure [4|. Another sharp decay in the free 
energy is found at around t = 25000. A closer investigation of the density field 
at t = 20000 and t = 28000 (Figure |4| shows that at t — 20000 the density field 
has a defect in the crystal lattice pattern which is relaxed by the removal of a few 
peaks (see t = 28000). The rapid decay in energy near t = 25000 corresponds to 
this process. This demonstrates the ability of the numerical method to capture 
the rich behavior of the hydrodynamic DFT model. 

An analogous simulation for the PFC model is run with e = 0.025, p = 
—0.07\/3 + \ and Cl = [0,32) x [0,32). The simulation is performed on a 
128 x 128 grid with s = 0.01. The evolution of the density field during the 
freezing transition is shown in Figure [6] The snapshots of the density with 
the flow field superimposed on it is shown on a zoomed region along side the 
density field on the full domain in Figure [6] The time evolution of the energy 
and its components are shown in Figure [7| It is once again observed that the 
total energy and the free energy are decreasing functions of time while the kinetic 
energy is non-monotone. Consistent with the nonlocal model, the kinetic energy 
again shows sharp increases in the regions corresponding to the rapid decay of 
the free energy and total energy. 
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6.3. Simulation of the Effect of Flow on Phase Transition in PFC Model 


In order to explore a practical application we study the effect of flow on 
a nanocrystal suspended in its liquid phase in a driven channel. The system 
considered is a channel represented by the domain Q = [0, L ) x [0, L), with L = 
25 6h where h = The boundary conditions are periodic in the horizontal 

direction and no slip moving wall boundary in the vertical direction: u^ n = 
{u W aii,n, 0 ) T and u i)0 = (u wa ii, 0 , 0 ) T at the boundary, where u wa ii,. corresponds 
to the horizontal velocity of the walls at the top (n) and bottom (0). No-flux 
boundary conditions for the density are imposed at the walls: pi, n = p^ n -1 and 


Pi, 0 — Pi, 1 • 

The crystal sample for the test is prepared by placing a crystal that is in 
equilibrium with the surrounding liquid and annealing it. The parameter e = 
0.4 is chosen where the co-existence region corresponds to pi = 0.72958 to 
p s = 0.90073. The solid phase under the one mode approximation takes the 


form 


p(x, y) := p s + A\ cos ( ffx + 


cos | -y 


- 2 cos (y) 


(76) 


where 



4\/3 

IV 


'I5e- 12 



(77) 


The initial sample is prepared by placing a solid phase given by Eq. (76) in 
a hexagonal region whose diagonal has the dimension 13a centered around the 
point (L/2, L/2). This hexagonal solid is surrounded by a liquid phase corre¬ 
sponding to homogeneous density pi. Here a represents the size of an atom 
a = An/y/3 which is the length of one period in x of the solid phase, one mode 
solution (Eq. (J76|). Thus the hexagonal nucleate shown in Figure [8] is 13 atoms 
wide along the main diagonal. The solid is expected to be thermodynamic equi¬ 
librium with the surrounding liquid. This initial condition is annealed until time 
reaches t = 20000, with wall velocity u wa n^/ n = 0 (see final density field in Fig¬ 
ure M) , where all components of the energy have equilibrated to a tolerance of 
10-“ All simulations in this test are performed at a step size of s = 0.02. The 
equilibrium field is not much different from the initial condition due to the spe¬ 
cific choice of p. During the annealing, however, the interface between the solid 
and the liquid becomes diffuse as expected. Thus the equilibrium configuration 
consists of a stationary crystal surrounded by a quiescent liquid. 

In order to study the effect of flow on the equilibrium crystal we now drive 
the system by moving the walls and setting up a parallel flow in the liquid. It 
is worth noting that the case where the wall moves with a non-zero velocity 
corresponds to a driven system. In this case the driving force from the wall will 
allow the system to be driven away from equilibrium. The system is no longer 
energy stable in the sense that the system will no longer follow the gradients of 
the energy to reach thermodynamic equilibrium. The results of the evolution at 
two different wall speeds (u wa n^/ n = ±0.1 and u wa u , 0 / n = ±0.5) are shown in 
Figures. |9| and [TT] respectively. The evolution in the low shear case shows that 
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the crystal does not change shape or size but merely rotates as the fluid flow 
is set up. A close-up of the velocity field that shows the rotation is presented 
in Figure [lOj However in the high shear case the crystal first begins to rotate 
and then change shape. As time progresses the crystal nucleate shrinks in size 
and settles into a steady smaller size. A closer investigation of the long time 
evolution (not shown) does indicate that the crystal stops shrinking beyond a 
critical size and continues to rotate in the fluid. This indicates that a system 
can be driven by flows to settle into a different phase co-existence than the one 
predicted by the stationary equilibrium phase diagram. This mechanism can be 
of critical importance in the growth of nanocrystals from a solution or liquid 
phase. 

7. Conclusion 

In this paper, we developed energy-stable, fully-discrete numerical methods 
for hydro dynamic-CD FT models in which compressible flow governed by the 
isothermal Navier-Stokes equations is driven by a free energy gradient. An ef¬ 
ficient nonlinear multigrid method was proposed and implemented to solve the 
implicit scheme. Although we did not demonstrate that the implicit discretiza¬ 
tion is uniquely solvable for any choice of time and space steps, the multigrid 
method was always able to solve the system without difficulty for a wide range 
of temporal and spatial grid sizes. 

Numerical simulations of both local (PFC) and nonlocal (CDFT) models 
were presented that demonstrate that the schemes are first order accurate in 
time and second order accurate in space. The energy stability and the ability of 
the methods to capture solid liquid phase transition in the model was verified 
numerically. Simulations illustrating the ability of the methods to capture the 
effect of flow on freezing transition and non-equilibrium properties of solid- 
liquid phase co-existence were presented. These simulations demonstrated the 
predictive ability of the model and the robustness of the numerical method in 
capturing physically relevant solutions. 

While first-order in time energy-stable methods, such as that presented here, 
are known to be significantly dissipative, we view this work as a necessary first 
step towards the development of second-order accurate (or higher) energy stable 
methods that would be significantly less dissipative. Nevertheless, the results 
presented here are accurate. The development of higher order accurate methods, 
as well as a more thorough investigation of the effect of flow using this model, 
and a Stokes flow counterpart, will be undertaken in future work. 

Acknowledgements. The authors gratefully acknowledge partial support from 
NSF Grants nsf-che 1035218, nsf-dmr 1105409, and nsf-dms 1217273. 
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Appendix A. Proof of Energy Dissipation in the General Hydrody¬ 
namic Model 

In this appendix we present a proof of Eq. <§• Note that by taking the time 
derivative of Eq. ([ 2 ]) we have 

= \ f dtp I u | 2 dr + [ pu d t u dr + [ ^-■ d t p dr. (A.l) 
at z Jq Jq Jq op 

Next we note that the hydrodynamic equations Eq. (|T]) can be written in prim¬ 
itive variable form as 


dtp + V • (pu) = 0, 


p ( d t u + u • Vu) = —pV ( — 


- qV 2 u. 


(A.2) 

(A.3) 


Using Eqs. (A.2) and (A.3) to eliminate the time derivative in the first and 
second integral we have 



(A.4) 

In the above calculation we have used integration by parts in steps 1 and step 
2. Further the cancelation in step 2 is achieved by the use of the continuity 
equation. 


Appendix B. Derivation of PFC model 

In order to derive an appropriate phase field crystal approximation we follow 
the work of Van Teeffelen et al [25] and Jaatinen et al [26]. First consider the 
CD FT free energy given by 


where 


and 


FcDFt[p\ = Fid[p\ + Fex [p ], 

(B.l) 

Fid[p\ := [ p(ln(p) - l)dr, 

(B.2) 

Jn 


Fex[p\ ■= J q P( J *P) dlc - 

(B.3) 
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The above free energy is approximated in two steps. 

Appendix B. 1 . Approximation of the Excess free energy 

Taking advantage of the fact that the dynamics is conserved i.e d t pdr = 0 
we can easily rewrite the excess free energy as 

Fex[p] ■= T^Pref jj,P~ Pref ) (J * (p ~ Pref )) dr, (B.4) 

where p re f is scalar value representing a homogeneous reference density field. 
Next we take advantage of the radial symmetry of the convolution kernel to 
expand it as a Taylor series in Fourier space about the zero mode: 


J(k) = C Q + C 2 k 2 + C 4 k 4 + ..., (B.5) 

where k is the Fourier variable. Now inserting the truncated form of the series 
into the free energy we have 

Fex\p] ■= \p 2 ref f 4>{C a - C 2 V 2 + C 4 V 4 )^r, (B.6) 

z Jn 

where 0 := p pref . Note that we have used the the real space representation 

Pref 

of the powers of the Fourier variables in terms of the gradients. Also note that 
this expression for the excess free energy is no longer non-local. The excess free 
energy now only depends on the local gradients of the density field. 


Appendix B.2. Approximation of the ideal gas part of free energy 

In a similar manner we will approximate the ideal gas part of the free energy 
as follows 


Fid\p] = Pref (4>+ l)[ln(0 + 1) + In {pref) - l]dr, 
Jn 


Pref 


f 

Jn L 


2 6 + 12 +--- + MPref)(J>+l) 


dr, 


(B.7) 


where we have Taylor expanded the logarithm term about 0 = 1. Note that if 
p re f is chosen to be close to average density there is reasonable expectation for 0 
to be small. Finally we note that the terms corresponding to In (p re j)(0 + 1 ) — 1 
do not contribute to the dynamics as the dynamical equations only depend 
on WSpElp]. This term does not alter the equilibrium as long as the mass 
conservation is enforced. Thus we drop these term and truncate the expansion 
to obtain the approximate expression 


id [p] 




dr. 


(B.8) 
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Appendix B.3. Phase Field Crystal Model 
Combining the two expressions (Eq 


B .6 


and 


B. 8 ) for the approximate free 


energies and re-arranging the terms we obtain the following expression 


J~\p\ Pref 


J 

Un 


' ~6 + u) dr+ 2 


2 + KqI 


V 2 ) 2 ]<^r 


(B.9) 


where the free constants and C 4 have been replaced without loss of 

generality with new constants A ,q Q and (. The above expressions can fur¬ 
ther be simplified by noting that the dynamical equations are conservative 
(i.e., dt f Q <pdx = 0 ) and that the dynamics only depends on the first vari¬ 
ational derivative of the free energy. Thus the terms which are multiples of 
<p 2 dr do not contribute to the dynamics and can be added and subtracted 
freely to complete powers. Now we introduce the change of variables r' = q Q r, 
ip = \J Iptr {fp ~ |) and e = ^4 (^ — C) and dropping terms that are linear in 
(p without loss of generality we have 





-e + (1 + V 2 ) 2 ]V> + 


dr'. 


(B.10) 


Note that the above expression is the standard PFC free energy [24] or the stan¬ 
dard Swift-Hohenberg free energy 031 . Through appropriate non-dimensionalization 
we can choose the constants q Q = A = p re f = 1 ([26]). For this choice of q Q ,A 
and p re f the phase and density fields are related by: 



and the approximate free energy can be written as 


(B.ll) 


Fpfc [f>\ 


= 3 [ % [—e + (1 + V 2 ) 2 ] ip + tpdr, 



|Vp| 2 + \(Ap) 2 


dr. 

(B.12) 


where a := 1 — e. This is the expression used this work. Also in the above 
equation the ' is dropped from the scaled r'. 


Appendix C. Proof of Convex Splitting Estimate 

Before we prove the Theorem |3.4| we need to obtain the following estimate. 

Theorem Appendix C.l. Consider a free energy functional F[p\ : H C 
L 2 er {Cl) —)> R, where H C L 2 er {Cl) is a Hilbert space of sufficiently regular 
periodic functions p : £4 —>• R with Q = [0, L x ) x [0, L y ). Then 

0-^)2 (C.l) 
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if 


^2 Jr \P + e^] > 0 


Vp, veH.ee 


(C.2) 


where (•,-)2 ^ the usual L 2 inner product and S p represents the variational 
derivative with respect to p. 

Proof. First we define a function G v , p {e) : M M for a fixed p.v e H as 
G v , p (e) := F[p + ev]. Now let us consider the case where Eq. (C.2) is true, i.e, 


G^ jV (e) is a convex function of e for all ip.v e H. From the convexity of G^ iV (e) 
we have 

Gv^(S)-G v ^( 0) = F[i/j + 5v]-F[i/>] >8Vim^-F[ip + ev] W>, v e H, 5 e M. 

(C.3) 

This gives us the inequality 

F[< i/> + Sv] - F[^\ > (S^F, ip + Sv — ^) 2 v e H. (C.4) 

where the defining relation for the variational derivative (in the Frechet sense) 


lim ~d e F ^ + = > 


(C.5) 


is used. Now identifying ip + Sv = <p in Eq. (C.4) we have the relation in Eq 


(C.l), this completes the proof of the theorem. 


It is worth noting at this point that Eq. (C.3) must hold for all S e M. This 
requires Eq. (20) to hold and not just the limit as given in E q. ([2l) ). Since 
we prove Theorem |3.4| using Theorem | Appendix C. 1 [ Theorem |3.4| requires a 
proper convex splitting (defined in Definition. 3.3) rather than a simple convex 
splitting. 


Now we are ready to prove Theorem 3.4 


Proof. (Theorem 3.4) 


Using Theorem | Appendix C.l| we have 

F c [ip\ - F c [(p\ > (S^Fcl^.ip - (p) 2 

and 

Fe[(p\ ~F e [lp] > (S^F e [lp].(p-lp)2 . 

Adding the two inequalities above we obtain the result. 


(C.6) 

(C.7) 


Appendix D. Simplification of Energy Dissipation rate 

In this appendix we prove the identity 

[ u fe+ 5 • V-D fc+ 5dr = -- [ V k+ i :V k+ ^dr. 

Jn 2 J n 


(D.l) 
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Integrating by parts and using the symmetry of V we have 


J u fe+ 5 • V-X> fe+ 5dr, =-\J (Vu fc+ 5 + (Vu fe+ 5) r ) : V k+ ^dr 

= -\J (Vu fc+ 5 + (Vu fc+ 5) r ) : (Vu t+ 5 + (Vu fc+ 5) T ) 
+ 1 J (Vu*+5 + (Vu fc+ y) : (V • u fe+ h) 


dr 


dr 

(D.2) 

where we have used the definition of V in Eq. (36) to obtain the second equality. 

It is easy to see that 

J (Vu fe+ 3 + (Vu‘+5) T j : (V • u fe+ 5l) dr = J (V • u fc+ 3l) : (V • u fc+ h) dr. 

(D.3) 

Now adding and subtracting - J (Vu fe+3 + (Vu t+5 ) T j : • u fc+ 2 l) dr to 

the right hand side of Eq. (D.2) and using the expression in Eq. flp^i ) we have 

f u fe+ 5 • V • V k+1 ?dr, = -I f (Vu fe+ 5 + (Vu‘+ 5) t ) : (Vu‘+5 + (Vu fe+ 5) r ) dr 
J O J o 

+ j (Vu t+ 5 + (Vu‘+ 5) t ) : (V • u fe+ h) dr 
f (V • u fe+ h) : (V • u fc+ h) dr, 

j (Vu fc+ 5 + (Vu‘ + yj : (Vu t+ 5 + (Vu H ^) T j dr 

y O 


1 

2 

1 

2 

-2 


j (Vu fc+ 5 + (Vu fc ^) T j : (V • U fc+ h) 

’J o _ 

V-u fe+ h) : (v-u fc+ h)dr 


dr 


= -- f : V k+ ^dr 

2 Jn 


This completes the proof. 


(D.4) 


Appendix E. Discrete Function Spaces, Operators and Inner Prod¬ 
ucts 

Appendix E.l. Discrete Functions Spaces 

In this section we define the discretization of the domain Q := (0 ,L X ) x 
(0, L y ). We use a nodal discretization. The development of the discretization is 
the similar to the one used in (Hang) where a staggered grid representation 
was used. 

Let h > 0 be the grid spacing such that L x = m • h and L y = n • h where 
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ra, n are positive integers. Now we define the sets 
Cm = {i • h | i = 1,... ,ra} 

Cm = {i • ft | i = 0,..., m + 1} (E.l) 

E rn — — 2 ? • • • j ^ T 2 } 

The set C m x C n partitions into a uniform rectangular grid with cells of size 
h x h. The set {Cfh x Cn}\{C m xC n } contains the ghost points outside the 
boundary of which are mapped back into in the case of periodic boundary 
conditions. 

Now we are ready to define the functions spaces as follows 

Cmxn = {0 : Cm x C n —» M} , Cjnxn = {0 : CW x CW > (E.2) 

Cfnxn = {0 • Cm X C n —>• M} , C mX n = {0 • Cm x CW , (E.3) 

Cxn = { W :4xC n 4M},Cxn = {^^x^4M} 5 (E.4) 

^rnxn = { u : x CW > Cxn = : C m x • (E.5) 

The spaces C mxn , Cmxn, Cmxn and C^xn contain the grid functions. A grid 
function are identified as faj := cj)(xi,yj ), where xi = i • h and yj = j ■ h and 
i, j are integers. The east-west edge-centered functions, in spaces and 

^rnxn are identified as Ui+ij •= u(x i+ i,yj). In a similar manner the north- 
south edge centered functions, in spaces £^f xn and £^ xn are identified as 
v i,j +i : = 

By defining a set of difference operators and equipping the space with a 
set of discrete norms one can write down self consistent, discrete gradient and 
Laplacian operators. 

We start by defining the relevant difference and average operators on the 
space. 

1. The averaging and difference operators A X ,D X : Cfhxn —> £mxn an d 
Ay,D y : C mX n —* £mxn are defined as 


•r.= !(/,.„ • /;../■)• Dxfi .! ( =* '!. 


Ayfij4* | 2 


(E.fi) 

2(/ij+i + Dyfi,j+± = j l (fi,j +1 — Aj)> j = 0 n 


The edge to center difference operators are defined as : £' 
and : £^ xn —^ C mxn are defined as 


’ew _v 

mxn f u mxn 


dJij = p/, 


)> dyfi.j ~ 


1 “ //.) i )‘ 


* = 1,..., m 

j = 

(E.8) 


2. The discrete gradient operator V/, : Cmxn xC mxS —5- C m xn x C mx « is 
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given by 


n+$,j + D x < t > i-\,jT D y ( t > i,j +1 + 


) T . 


i = 1 ,..., m 


^h^ij — 2 

3. The discrete divergence operator is defined in a similar manner as 


(E.9) 


V/i • u i:j — ^ [D x u i+ i j + D x u i _i j + D 


y v i,3+\+ D y v i,3-\)' 

i = 1,..., m 


(E.10) 




where we have used the definition u^- := (uij,Vij) T . 


4. The discrete 2 D Laplacian A/,, : Cmxn —> C m xn is defined as 

i = 1,..., m 


^h^ij - d X {D x (j))ij T dy(Dy(f)^ij 

= ^2 (0*+l>J + 0i,J+l + fa-1, j + fa-1,j - 4 fa,j) j = l n 

(E.ll) 


We define the following weighted inner products on the function spaces 

m n 

(0 II VO = h ^ ^ ^ ^ ■> ^ Cmxn U Cffixfi U Cmxn U Cmxn (E.12) 

i—1 j =1 

with u = (ui,U 2 ) t and v = (fi, the inner products 

(ll || v) = ( U\ || V\) + ( U 2 || ^ 2)5 ^ 1 ? ^ 2 , ^ 1 5 ^2 £ C m xn U Cmxn U C m xn U Cmxn 


i=l j=1 
m n 

[/llfi'lns = ^ E E (/ij+ift,j+i + J 

i=l J = 1 

and the one dimensional inner products are defined as 

n 

(f*j+h I 5*,j+i) = h E fi,j+k g i,j+ 5’ e ^ 




(E.13) 

•j g < - ,].j) ' 

> fi9£ £ mxn 

(E.14) 

5 ^ 0 - 5 ) ’ 

f,9& £mxn 

(E.15) 


z=l 


ns 

mxn 


(/z+i,* I d 22 fi+ ~ ,j9i+ 1,j 5 f 1 9 ^ 




! ew 

mxn 


0 *,j I ^*,j) = h J2 ^ec v 


(E.16) 


Z=1 


Oz,* | ^i,*) = h E ^ e Cmxn 


z=l 
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Using the definitions given in this Appendix and in El, we obtain the fol¬ 
lowing summation-by-parts formulae: 


Proposition Appendix E.l. (summation-by-parts) If (f G Crn Xn UCrnxn and 
f e £™ xn then 


[Ar<A||/] ew 


- {<t>\\d x f) 


i^A x (j) l* /i,*) + (A^, 


'ra+|,* 



and if<pe C m xn U C SX n and f € £“ xn i/ien 


(E.17) 


[A/ 0 II/L 


OAK/) 


(aA,^ /*,§) + (A/<a. 


T*r,71-|- 2 



(E.18) 


Proposition Appendix E.2. (divergence theorem )<fi G Cmxn and u = 

(it, 1 ?) G Cmxn X Cmxn - Then 

(V^ || u) = — (0 1 V/i • u) 

2 (00,* I ^1,*) A 2 I ^ra+1,*) ^(01,* I ^0,*) + ~(^m+1,* | / a m ,'k) 

— 2^*’° I ^*> 1 ) 2 I ^*, n +i) — 2 (0*4 I ^*, 0 ) d~ 2 ^*’ n+1 I ^*, n )‘ 

(E.19) 


Proposition Appendix E.3. (discrete Green’s first identity) Let q i>, if G Cmxn- 
Then 


[D x (j>\\D x ip\ ev/ + [A^HA/V’hs = -(^IIA^) 

— ^^01,* Dx^/ji^ + ^Aaj0 m _|_i ^ -DxV’m+U*) 

— + ^A/ 0 * 5 n+J ^ 2 / 0 *,n+§) • 

(E.20) 

Proposition Appendix E.4. (discrete Green’s second identity) Le£ </>, ijj G 
Cmxn- Then 

(011 A/i^) = (A/, 011 ^) 

+ (Ac0 m +|,* _ (^E0ra+J,* ^VVn+U*) 

— ^ 2 C 0 i,* H" ^EU 0 i ? * 


+ ^A y (^ n+ i ^3/0*,n+5) (^2/0*,n+§ A/0*,n+|) 

— ^Ay (/>*,! Dy'ljj+ i'j + ^-^3/0*, I Ayfj^ l^ . 


(E.21) 


28 


















Appendix E.2. Boundary Conditions 

In this paper we work with periodic functions. We will consider functions 
0 £ Cnxm which have the form 

0m+ l,j = 01, ji 00, j = 0m, j 5 j = lj • • • , n, 22) 

02,72+1 = 02,1, 02,0 = 02,72, i = 0, . . . , 777. 

For such functions the center to edge averages and differences are also periodic. 
Further at this point we note that the results proven for periodic functions in 
the following sections will also hold with slight modifications for homogeneous 
Neumann boundary conditions: 


0m+l,j — 0m, j , 00, j — 01, j 5 j ~ 1? • • • 5 23 ) 

02,72+1 02,725 02,0 02,15 ^ O 5 • • • 5 

and homogeneous Dirichlet (no-slip) boundary conditions for velocity, e.g. u = 
(^1,^2) = (0,0): 


^lc,m+l,j — 'U j k,m,ji ^k,0,j — ^/c,l,j 5 j — 1 , . . . ,77^ k — 1 , 2 

^k, 2,72+1 = ^/c,2,72 5 ^/c, 2,0 = R'/c + l, 7 = 0,... , TYl] k = 1, 2 


(E.24) 


In the periodic, Neumann boundary and no-slip condition cases, the bound¬ 
ary terms from the integration by parts vanish. 


Appendix E.3. Discrete Norms and Their Properties 

We now define norms for the cell centered periodic functions. If 0 G C mxn 
and periodic then 

|| <j> || 2 := +(++), (E.25) 

|| cP 11 4 5~ s/W II !)• (E.26) 


Appendix F. Proof of Theorem 5.1 


In this appendix we present the details of the proof of Theorem 5.1 


Proof. Summing Eq. (48) over all cell centers and using the discrete Diver¬ 
gence theorem, we immediately have mass conservation (p /c+1 || l) = (p fe || l). 
Choosing 0 = p k+1 and -0 = p k in Eq. (53) we have the following estimate on 
the free energy 


+/ +1 ] - +/] < (++[/ +1 ] - 5 p F e [p k },p k+1 - p k ). 


(F.l) 


This estimate immediately translates to an estimate for the discrete total energy 
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as follows : 


£\p k+1 , u fe+1 ] - £[p k , u fe ] = \{p k+1 n k+1 || u fc+1 ) - X - (p k u k || u fe ) 

+T[p k+1 ]-T[p k } 

< \{p k+1 n k+1 || u fe+1 ) - 1 (p k u k || u fc ) 

+(p k+1 -p k || p k+l ), 


where we have used Eq.(40) in in step 2 along with the definition of p k+1 (Eq. 


34). 


Multiplying Eq. (49) by u fc+ 2 and summing over cell centers we have 

\{p k+1 u fc+1 II u fc+1 ) - l -{p k n k || u fe ) -\{{p k+1 -P k ) III u fc+1 I 2 ) 

= -s l -{p k u fe +3 || v h I u fe+1 I 2 ) (F . 3) 

-s(y u fe+ 51| v„/ +i ) 

+s 7 (u fe+ 5 || A h u k+ ^), 


Now simplifying this in a manner analogous to the proof of Theorem 4.1 
have 


1 


^{p k+1 u /c+1 || u /c+1 ) - (p k u 


1 


k„k II u k\ 


+(p /c+1 - p k ,p, k+1 ) = +S7(u fe+ a II A/,u fc+ 2 


— sy ([Th^ll Ar^] ew + [ D y v\\D y v ] ns ) , 

(F.4) 

where u : u are components of u. Here we have used the discrete Green’s identity 
(see Proposition Appendix E.3) in the last step. Finally using the relation Eq. 
« and Eq. (F.2) we have 


£[p k+1 ,u k+1 }-£[p k ,u k ] < \{p k+l u fc+1 || u fc+1 ) 

-\{p k n k \\ u k ) + (p k+1 - p k |] p k+1 ) (F.5) 

= -sj ([D x u\\D x u] ew + [D y v\\D y v} ns ) 

< 0 . 

Hence we have shown that the total energy of the system is non-increasing, 
regardless of the time step s > 0 i.e £[p kJrl , u /c+1 ] < £u k ]. 


Appendix G. Nonlinear Multigrid Solvers 

In this Appendix we present the details of the nonlinear multigrid algorithm 
used for time stepping the implicit numerical scheme presented in the previous 
sections. 
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Appendix G.l. Multigrid method for the convex splitting scheme for the CDFT 
model 


In order to time step the implicit numerical scheme presented in Eqs. (48) 


- (51) with the CDFT free energy Eq. ( [56] ) ? Eq. (58) and Eq. (59) we need 
to solve for p /c+1 , u k+1 , v k+1 G Cfhxn such that = (u k + 1 ,v^ 1 ) T 

p k ,u 


given 

' k ,u k ,v k G Cfhxn such that u k = (u k ,v k ) T , that satisfy the following equations 

4“ (Au Ai - rUX-l, + A+,AA - A-.Ai) 


(G.l) 




k k k-\- h 

-Y h (i AA i 2 -1 “A i 2 ) 

P (| „.fc+l 12 


pi 


v: 


i+lj 


~ V, 


fc +1 

*-l,j 


2 P«(^e*l) / fc+1 _ fc+1 \ 

V^i+lj' Pi-1,j) 


2 h 


+ ^( J *Pi + U- J *Pi-l,j) 

+ ‘FAh 12 (^ +1J - ^_ u ) + 7A /1 u: 


(G.2) 


fc+1 

+7 


and 


Pij (v 


y 1 tj 


v k ■) = 

U ljJ 


where 


c d, 


k k fc + h 

PijUijUij 

4 (K fe+1 


, 2 - I u k+1 
d+il I 


(I Aii i 2 -1 A-i i) 

“ if ( ln (Aii) - ’“A-i» 

2 Pij+e*l) / fc+1 fc+1 \ 

2h Vfc’+i+i Pi,j-i) 

+ §d*4 J+ l-^*4+l) 


2 h 


2h A« 




(4j+i - 4i-i) + ++4 


4m + u i,m) • 


fc+i 


(G.3) 

(G.4) 
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The above non-linear problem can be written in terms of a non-linear oper¬ 
ator N and source S such that N(g) = S. Let g = (p, u : v) T then the 3 xmxn 
nonlinear operator N(g) = (N^\g)^ N^ 2 \g)^ (g)) T can be defined as 


^ij\ g) Pij + ^ (Pi+l,j u i+lJ Pi-1,j u i-l,j + Pi,j+l v i,j+l Pi,j-l v i,j-l) ? 

(G.5) 


Ng } (g) -.= p%u.. 




2 Pij u ij Vi i 


~4h ^ Ui+hj ^ “ I Ui ~ ld ^ 
(I v i+hi I 2 - I Vi-U I 2 ) 

2 rf S (Je*l) . , 

2h (Pi+i,j Pi-i, j) 

+ 7^A.hUij 


N « } (g) : = Pii v ij ~ s 


2 

p k 




u i,j~ 1 


_ 4 | (I Vi ’ i+1 I 2 “ I ^ 


■^T- (MPij+i) -ln(Pi,j-i)) 

2/4-(Je*l) 

2h vPhj-i Pij-i) 

~^2^h v ij 


(G. 6 ) 


(G.7) 


where uo^. is as defined in Eq. (G.4). The 3 xmxn source S = (S^\ S^ 2 \ S^) 1 


is defined as 


o(!) — n k _ ( k k _ k k , k k _ fc k \ 

n ij Pij ^ \Pi+l,j a i+l,j Pi-1,j a i-l,j ^ Pi,j+l u i,j+l Pi,j-l u i,j-l) i 


(G. 8 ) 


s£?> := 


• — U — ?/ — 
rij a lJ 




rij ^ij u ij 


+^( J *Pi + i,j-J*ptu) 

2 4(J< 


+ 


1) 


2 h 


(Pi+ 1,3 


Pi-1, j) + 9 A h Uij 


(G.9) 
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■■= PiA +s 


-n k .M k .v k . 
) PlJ WIJ U IJ 

h 

P\ 


+ %{J*pl j+ i-J*fo- 1) 


2 h 

2,44*1) 


(G.10) 


2h (Pi,j +1 - Pij-i) + 2 


Given the Gth guess (//, ?/, i/) for (p /c+1 , ^ /c+1 , ?; /c+1 ) we can get the (£+1)- 
th guess for the (p kJrl , u kJrl , v k+1 ) using the following smoothing scheme. The 
smoothing scheme takes the form of a 3 x 3 linear system which is in turn solved 
exactly using Cramer’s rule. 


J +1 = Q (1) 4- A (n k 7 / - n k 7/+ 1 4- n k v £ - n k 7/+ 1 \ 

Pij ^ij + 4^ \Pi+l,j U i+l,j Pi-l,j U i-l,j ' Pi,j+l V i,j+l Pi,j-l V iJ-l) 


(G.ll) 


0 H. + ] u *+ 1 + -n k -u k V +1 

Pw ' 2h, 2 / o ' 


= S 


( 2 ) 




--(I 

/I h V 1 


4h Vl <+1 ’* 


- u, 


W 


i-fhj 


~ V\ 


i 2 ) 


ji 2 ) 


p K 

~ (Hpl+tj) -Hpi-l,j)) 

2Pij(Je * 1) , e p + 1 . 

4+1,J Pi-l.jJ 

-o&(‘ 


7 2/1 

2/? 


■/+ 1 


u, 


hi+i 


<th) 

(G.12) 


4 + ) ^+! _ - 0 fe . w fc +/ +1 


= S 


(3) 






- U, 


r+i 


4+i |2\ 


-1 I 2 ) 


( in (4i+i) - in (4t-i)) 

^Pij{Je * 1) / p _ p+ 1 X 

2h V^,i+i Pi,j-i) 

+ 2 ^2 (4+1 j + 4-h + 4/ + 1 + <i-i) 

(G.13) 
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Appendix G.2. Multigrid method for the convex splitting scheme for the PFC 
model 


In order to time step the implicit numerical scheme presented in Eqs. (48j) 
(51) with the PFC free energy in Eqs. (63)-(65)we need to solve for p /c+1 , 


v k+1 G 


Crhxn such that u /c+1 = (u k+1 , v k+1 )' J given p k ,u k ,v k G Cfhxn such that u k 
(■ u k ,v k ) T that satisfy the following equations 


p k+1 

rij 


p k . = 


T h C 


k+b k k +h , k k +h k 

- Pi,j- 


Pi+ l,5'“t+l,5 


j-i v i,A) > 

(G.14) 


I’M)' 1 ~ u ii ) = 


it 


fe+i 


k k &+ ^ 

-O-UJ-V- • 2 

rij lj IJ 

-A (\ v k +1 |2 _ 

4£ ^ i+1 ’ J ' 1 

-fj (I 4+ + ft I 2 - I »f-w 12 

-§Wi‘ J '-Kfi', J ) + 74»»: 


i—1,5 I ) 

) 


fe+i 


(G.15) 


4(44-4) 


J7 Jo /c+ A 


4 

4ft 

4 


7, /c + 1 

^4 + 1 


l i,3 -1 




/c+l 


hi+i 

fc+i 


*,5 


41 2 ) 






(G.16) 


and 


fc+i = £ fe+i _ £ , a 0 fc+i _ £ 

3 2/ + 2 


+ 2Afjpk. -)- A/jk: 


/c+1 


(G.17) 


where 


4 +1 = A ^4 + t 


(G.18) 


4 = 2/j (4+i,j - 4- u - 4,j+i+4i-i) • 


(G.19) 


The above non-linear problem can be written in terms of a non-linear opera¬ 
tor N and source S such that N = S. Let g = (p, iq i>, /q k) t then the 5 x m x n 
nonlinear operator N(g) = (TV^^g), iV^(g), iV^(g), 7V^(g), Af( 5 )(g)) T can 
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be defined as 


-ivt(I) / \ ° f k k k 

N ij (g) : = Pij + jr (Pi+i,jUi+i,j ~ Pi-i,jUi- hj + Pi , j+ 1 v iJ+ 1 - p, 


4 h 


*,3- 


N ij (g) : = Pij^J 


Ng'(g) : = Mr - 5 (ft, 

r(5) 


—s 


\pij^ijVij 

^ (I U i+hj I 2 - I u i-l,3 I 2 ) 

-^(1 v m,j I 2 - I I 2 ) 

+ 2^ (Mi+iJ _ Pi-1,j) 

~^~2 ^h u ij 


N ij (g) : = 4%' 


Pij^ijUij 


r_ 

4 


u i,j +1 I I u i,j-l I ) 

v<j+i I 2 - I I 2 ) 


+ 2^ (Pi,j -1 _ Pi,j- 1) 

+ 2 ^ hVi 3 


k +1 


— a 


N^(g) := /c^-Afcp; 


4 +1 -£)- a *4 +1 ’ 


n /c+l 


where uj k , is as defined in Eq. (G.19) 


The 5 x m x n source S = (. SS^ 2 \ S^ 3 \ S^\ S^) T is defined as 
S S } := Pij ~ 4 (Pi+l,3 U i+l,3 ~ Pi-l,j U i-l,j + Pi,j+l V i,j+l ~ Pi,3- 1 1 


o(2) , = k k 
ij ’ rij ^ij 


Sg } :=44 


-^4<44 + ^ A *4 
^444 +| A >4 


sg } := 2A fc p& 

sg } := o 

Given the Gth guess (//. </. r f .//. f/) for (p k+1 , u k+1 , v k+1 , p k+1 
can get the (£ + l)-th guess for the (p k+1 , u k+1 , v k+1 , p k+1 , n k+1 ) 


1 v i,j 0 > 

(G.20) 


(G.21) 


(G.22) 

(G.23) 

(G.24) 

,,k \ 

i,j- 1/ ’ 

(G.25) 

(G.26) 

(G.27) 

(G.28) 

(G.29) 

, K k+1 ) we 
using the 
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following smoothing scheme. The smoothing scheme presented here takes the 
form of a 5 x 5 linear system which is in turn solved exactly using Cramer’s rule. 




— Cl 1 ! I_ ( n k _ n k ? /+l JL n k n/ _ n k «/+* I 

°ij + \Pi+l,j u i+l,j Pi-1,j u i-l,j Pi,j+l v i,j+l Pi,j-l u i,j-l) > 


(G.30) 


0 k. , ^L] u w , i 0 k k t +1 _ 
P 13 ' 2/^2 y ‘hj ' 


( 2 ) 






p_ 

4 i 


'c; 


i+X,j 


- V. 


/+1 
1 J 


T/ _ 7 /+l \ 
^ 2h v i + 1 >J Pi-i,j) 

7 , ..€+i 


2^2 ( U *+l,i “ r 


l i,j+1 


+ <t -0 

(G.31) 


+ Td ^+1 _ - J O fe .M fc .^+ 1 

j ' 2 / 1 2 / *•'/ 2 ■' u i:i 


P_ 

4 h 


u i,j+l 


S^+s 

+ 2 h 2 ( Vi+1 ’ j 


^i,j-\-l 


I 2 ) 


- V. 






— // +1 
Pi,j-l) 


4i- 


4,7+1 


(G.32) 


(a + (pC - 3/2) 2 ) f /.j 1 


+ ,/ +1 4. —k €+1 

^Pij ^ ^2 u 


a (4) 

2 

3 


4 


(a + (pC - 3/2) ) 

+ 7^ G+ij + + di+1 + d,j-i) > 


3\ 3 

2 ) + 2 


(G.33) 


ft!Pip + d/ 1 - S if + 7^2 (Pi+hj + j°f—+ Pi,j+1 + Pij-1 ) • (G.34) 


ft 2 


In Eq. (G.33) we have linearized the cubic term using the approximation 
( 0 ^+i ) 3 po 3 (^) 2 ^+i _ 2(<^/) 3 within the lexicographic index C 
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Appendix G.3. Damped Gauss-Seidel Iteration 

To further aid the convergence of the method and improve the stability of 
the Gauss-Seidel method we use a damped Gauss- Seidel method. Given p k and 
pk+i,t g Uess for p k+1 , let y r/ c + 1 J+bG'S' be the {£ + l)-th guess obtained 

from a smoothing scheme. The damped Gauss-Seidel method updates the new 
guess p /c + 1 ^+ 1 as 

p k+u+i = (1 _ Wdamp ) p k+U+1.GS + Wdamp pk+w (G.35) 

where Wdamp is a damping co-efficient set to be Wdamp = 0-5 in all simulations 
presented. 
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h f 

II || 2 

Rate of Convergence of p 
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4.8685 x 10“ 8 

2.3576 

32 

128 

32 

256 

1.0802 x 10 -8 

2.1722 
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256 

32 
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2.2465 x 10“ 9 

2.2655 


Table 1: The error and convergence rate for the density for the convergence test of the 
hydrodynamic CDFT model discussed in Section [6.1.1| 
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Rate of Convergence of u 
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1.0715 x 10“ 8 
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256 
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2.5396 x 10“ 9 

2.0770 


Table 2: The error and convergence rate for the horizontal component of the velocity for the 
convergence test of the hydrodynamic CDFT model discussed in Section [6. 1.1 1 
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Table 3: The error and convergence rate for the vertical component of the velocity for the 
convergence test of the hydrodynamic CDFT model discussed in Section [6. 1.1 1 
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Rate of Convergence of p 
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Table 4: The error and convergence rate for the density for the convergence test of the 
hydrodynamic PFC model discussed in Section [6.1.2| 
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Table 5: The error and convergence rate for the horizontal component of the velocity for the 
convergence test of the hydrodynamic PFC model discussed in Section [6.1.2| 


h c 

hf 

II e h v ^ || 2 

Rate of Convergence of v 

32 

16 

32 

32 

6.6965 x 10 -7 

— 

32 

32 

32 

64 

1.4285 x 10“ 7 

2.2289 

32 

64 

32 

128 

3.2625 x 10 -8 

2.1305 

32 

128 

32 

256 

7.7450 x 10 -9 

2.0747 

32 

256 

32 

512 

1.8833 x 10“ 9 

2.0400 


Table 6: The error and convergence rate for the vertical component of the velocity for the 
convergence test of the hydrodynamic PFC model discussed in Section [6X2] 
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Figure 1: The evolution of kinetic energy, Helmholtz free energy T[p] and total energy £[p, u], 
as labeled, for the simulation cor respon ding to the convergence test of the Hydrodynamic 
CDFT model described in Section [6.1.1| with (h = 512). 
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Figure 2: The evolution of kinetic energy, HelmholtzfFree energy T[p] and total energy £[p, u], 
as labeled, for the simulatio n corr esponding to the convergence test of the Hydrodynamic PFC 
model described in Section [6.1.2| with (h = 512). 
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Figure 3: The figure shows the disc rete Fourier image of the convolution kernel J used in the 
freezing simulation (see Section [6.2| >. The Fourier image is radially symmetric and represented 
in 2-D (right) and in 1-D (left) along the radial direction. 
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Figure 4: The time evolution of the density field with the velocity field superimposed. The 
figure shows the liquid to solid phase transition in the hydrodynamic CDFT model discussed 
in Section EJ 45 
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Figure 5: The evolution of kinetic energy, Helmholtz free energy T\p\ and total energy £[p, u], 
as labeled, for the simulation of freezing of the hydrodynamic CDFT model described in 
Section [6~2] and shown in Figure^] 46 
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Figure 6: The time evolution of the density field. The figure shows the liquid to solid phase 
transition in the hydrodynamic PFC model discussed in Section [6.2| 
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Figure 7: The evolution of kinetic energy, Helmholtz free energy T[p] and total energy £[p, u], 
as labeled, for the simulation of freezing of the hydrodynamic PFC model described in Section 
|6.2| and shown in Figure [6] 
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Figure 8: The initial and final configuration during the annealing process of a nano-crystal 
placed in a channel with d ensit y of solid and liquid phases chosen in the co-existence regime 
at equilibrium (see Section |6.3| >. 
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jure 9: Snapshots of the density field at different times are shown as a nanocrystal is sheared 
a driven channel with wall speed u wa u ^ n /o = ±0.1 as described in Section 16.31 The initial 
iditions correspond to a system at equilibrium with u = 0 and u wa u Q/ n = U (rig. [8]). The 
ire shows that the crustal tumbles in the channel while holding its shape fixed. The colors 
ailable on-line) correspond to a linear RGB scheme from 0 to 2.9 as shown in Figure [8] 
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Figure 10: The flow field superimposed on the density field over a zoomed-in region of the 
domain at times t = 3200.0, 3600.0 and t = 4000 from the simulation presented in Figure [9] 
The figure shows the velocity field does not correspond to that of a plane Couette flow but 
changes to accommodate the tumbling nanocrystal. The colors (available on-line) correspond 
to a linear RGB scheme from 0 to 2.9 as shown in Figure [8] 
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jure 11: Snapshots of the density field at different times are shown as a nanocrystal is 
:ared in a driven channel with wall speed u wa u , n / 0 = ±0.5 as described in Section|6.3| The 
bial conditions correspond to a system at equilibrium with u = 0 and u wa u f Q/ n = 0 (Fig. 

The figure shows that the crystal initially shrinks to a smaller size and then stabilizes in 
3 as it tumbles in the channel. The colors (available on-line) correspond to a linear RGB 
erne from 0 to 2.9 as shown in Figure [8] 
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